R语言元分析专题第四章：合并效应量

第四章合并效应量(Pooling Effect Sizes)属于每个元分析的核心部分，合并或汇总效应量是为了获得研究的整体效应大小。

本章讲述对已整理的数据文件，在R中实现效应量合并分析，而对单个研究效应量的提取可见第十五章内容。对研究总体效应量的估计存在两种方法：固定效应模型(Fixed-Effect Model)和随机效应模型(Random-Effects Model)。

两种模型的前提假设不同，固定效应模型假设所有的研究及其效应大小都来自单一的同质群体(Borensteinet al., 2011)。反之，异质群体选用随机效应模型。关于哪一种模型最适合哪一种情况仍存在着广泛的争论(Fleiss,1993)，在临床心理学和健康科学(Cuijpers, 2016)中建议只采用随机效应模型，以下将讲述R中两种模型代码的实现。

使用meta包(Schwarzer,2007)。

library(meta)

4.1 Fixed-Effects-Model

4.1.1Pre-calculated effect size data

为了计算整体效应量，需要对所有效应大小进行平均，但对于精确研究给予更高的权重。在这种情况下，更高的精度意味着研究有更大的样本量N，评估效应量时才会有更小的标准误差（Standard Error）。

在第3.1章，Excel电子表格有两种方法整理元分析数据：它可以存储为原始数据（包括每个研究的平均值、N和SD）；或者它只包含计算出的效应量和标准误差(SE)。

根据使用数据格式的不同，合并固定效应模型的函数也有所不同，但差别不大。对于预先计算好的效应量数据（如madata数据集），其中包含每个研究效应量和SE。

该数据中的效应量是基于连续的结果数据。由于已经计算了效应量，可以使用meta::metagen函数。对于这个函数，首先可以指定参数进行加载，所有这些参数都可以在meta包加载后通过在控制台中键入?metagen来查询。

下表是代码中最重要的参数解释：

固定效应模型的元分析代码，对这个分析的结果简单地命名为m。

m <- metagen(TE,
             seTE,
             data=madata,
             studlab=paste(Author),
             comb.fixed = TRUE,
             comb.random = FALSE,
             prediction=TRUE,
             sm="SMD")
m

输出结果：

固定效应模型元分析的结果包括：

• 每个研究的单个效应量和它们的权重；
• 纳入研究的总数(k)；
• 总体效应量（g= 0.48）及其置信区间和p值；
• 研究间的异质性测量，如tau²或I²，以及异质性的Q检验。

使用$命令，还可以直接查看各种输出。例如：

m$lower.I2

给出了I²的95%置信区间的下限。

## [1] 0.3787897

可以使用以下命令将元分析的结果保存为.txt文件到我们的工作目录中：

sink("results.txt")
print(m)
sink()

4.1.2 Raw effect size data

原始数据（数据已经按照在3.1.1章中描述的方式准备好）使用meta::metacont()函数来替代3.1.1函数。这些代码的结构看起来非常相似。

将使用数据集metacont，它包含了所有综合研究的原始数据。

使用meta::metacont函数和metacont数据集。现在给输出命名为m.raw。

m.raw <- metacont(Ne,
                  Me,
                  Se,
                  Nc,
                  Mc,
                  Sc,
                  data=metacont,
                  studlab=paste(Author),
                  comb.fixed = TRUE,
                  comb.random = FALSE,
                  prediction=TRUE,
                  sm="SMD")
m.raw

输出结果：

4.2 Random-Effects-Model

在此之前，讲述了如何使用metagen和metacont函数执行固定效应模型的元分析。

然而，固定效应模型仅适用于假定所有纳入元分析的研究都来自同一个总体（即各项研究的真实效应值相同），而在实际情况当中往往不是这样的。比如，各个研究的真实效应值可能会随着各项研究中的干预措施的不同，样本、研究方法的不同而更大或更小。在这种情况下，我们不能认为所有的研究都来自于同一个“总体”（我们没有理由假设所有研究的真实效应是相同的）。

同样地，一旦我们在使用固定效应模型进行元分析检验统计异质性，显示I²>0%的时候，就提示我们应该使用随机效应模型去进行元分析（I²：研究之间真实异质性占总变异的比重）。当然，在R当中，使用随机效应模型也并不繁琐。

4.2.1 Estimators for τ2 in the random-effects-model

操作上，在R中进行随机效应模型元分析与进行固定效应模型元分析并没有太大的不同。然而，我们选择τ²作为估计量。下面是在meta中实现的估计，在元分析代码中的使用method.tau变量。

所有这些估计推导τ²使用稍微不同的方法，导致不同的合并效应量估计值和置信区间。这些方法或多或少的偏倚，常常取决于许多研究k的背景、参数和在每一个研究中被试数量n。n因研究而异，τ²随之变化。

Veroniki和他的同事发表了一篇综述性论文，对当前的证据进行了很好的总结(Veroniki et al., 2016)。文章可公开访问：

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4950030/

在医学和心理学研究中，目前最常用的估计量是DerSimonian-Laird估计量(DerSimonian & Laird,1986)。这种广泛使用的部分原因可能是RevMan或综合元分析（旧版本）等程序只使用这种估计。

然而，模拟研究表明Maximum-Likelihood、Sidik-Jonkman和Empirical Bayes估计在估计研究间方差方面具有更好的表现(Sidik & Jonkman, 2007; Viechtbauer, 2005)。

Hartung-Knapp-Sidik-Jonkman方法

对DerSimonian-Laird方法批评是，当我们的估计合并效应量的方差var()，该方法容易产生假阳性(IntHout, Ioannidis, & Borm 2014)。特别是当研究的数量很少，且存在很大的异质性时(Hartung, 1999; Hartung & Knapp, 2001a, 2001b; Follmann &Proschan, 1999; Makambi, 2004)。

因此，Hartung-Knapp-Sidik-Jonkman(HKSJ)方法是一种更稳健的估计var()。已经证明，在许多情况下(IntHout, Ioannidis, & Borm, 2014)，该方法显著优于DerSimonian-Laird方法。HKSJ方法也可以很容易地应用在R中，而其它程序还没有这个选项。这是在R进行元分析的另一大好处。HKSJ方法通常导致更保守的结果，以更宽的置信区间表示。

4.2.2 Pre-calculated effect size data

与固定效应模型第4.1章相比，随机效应模型元分析只需要定义三个额外的参数。正如之前所描述的，如果我们想使用Knapp-Hartung-Sidik-Jonkman调整，必须告诉R研究间方差估计量(τ²)。

下面是我们必须在代码中定义的所有参数表，这些参数用于执行具有预先计算效应量的随机效应模型元分析：

将再次使用madata数据集进行元分析。使用Sidik-Jonkman估计(“SJ”)和HKSJ方法。进行此分析前，请确保在R中加载了meta和metafor。

library(meta)
library(metafor)

由于效应量数据是预先计算的，所以将使用meta::metagen()函数。

m.hksj <- metagen(TE,
                  seTE,
                  data = madata,
                  studlab = paste(Author),
                  comb.fixed = FALSE,
                  comb.random = TRUE,
                  method.tau = "SJ",
                  hakn = TRUE,
                  prediction = TRUE,
                  sm = "SMD")
m.hksj

随机效应模型输出结果：

输出结果显示，估计的效应量是g=0.59，95%置信区间从g=0.39延伸到0.80（四舍五入）。与在第4.1章的固定效应模型元分析中发现的效应（g=0.48）相比，这种效应更大。

将其与使用DerSimonian-Laird估计输出进行比较，并在设置hakn = FALSE进行比较。因为这个估计是默认的，所以这次不必定义method.tau。

m.dl <- metagen(TE,
                seTE,
                data=madata,
                studlab=paste(Author),
                comb.fixed = FALSE,
                comb.random = TRUE,
                hakn = FALSE,
                prediction=TRUE,
                sm="SMD")
m.dl

输出结果：

可以看到，使用这个估计的总体效应量估计值与前一个类似（g=0.57），但是置信区间更窄，因为我们没有使用HKSJ方法对其进行调整。

4.2.3 Raw effect size data

如果使用原始的效应量数据，例如存储在metacont数据集中的数据，再次使用meta::metacont()函数。参数保持不变，这次把元分析结果命名为m.hksj.raw。

m.hksj.raw <- metacont(Ne,
                       Me,
                       Se,
                       Nc,
                       Mc,
                       Sc,
                       data = metacont,
                       studlab = paste(Author),
                       comb.fixed = FALSE,
                       comb.random = TRUE,
                       method.tau = "SJ",
                       hakn = TRUE,
                       prediction = TRUE,
                       sm = "SMD")
m.hksj.raw

随机效应模型输出结果：

4.3 Binary outcomes

在某些情况下，研究人员需要使用二分类结果数据（binary outcome data）（例如：死亡/存活、抑郁障碍/无抑郁障碍)，而不是连续的结果数据。有两种常见的二分类结果数据：

事件率数据（Event rate data）。在这些数据中，只处理每个组中经历事件的人数，以及每个组中的总样本量，从这些数据中计算出的效应量包括比值比、相对风险度或风险差等。

发病率数据（Incidence rate data）。事件率数据通常不包含关于事件发生或不发生的任何时间跨度的信息。考虑到研究的随访时间通常有很大的不同，因此考虑事件发生的时间间隔也很有意义。在流行病学中，发病率通常用来表示在一个标准的时间段内发生了多少事件。相应的效应量为发病率比（IRR），IRR是将干预组的发病率与对照组的发病率进行比较。

对于已计算出的效应量。在这种情况下，可以像以前那样使用metagen函数（参见第4.1章和第4.2章），还可以参照一些其他说明。将在第4.3.3章中描述如何使用metagen函数来处理预先计算的二分类结果数据。

只有原始的结果数据。如果是这种情况，我们将使用meta::metabin()或meta::metainc()函数来代替。

4.3.1 Event rate data

对于事件率数据的元分析，使用数据集binarydata。

load("_data/binarydata.RData")
str(binarydata)

进行随机效应模型元分析，汇总度量值是风险比（RR）。

m.bin <- metabin(Ee,
                 Ne,
                 Ec,
                 Nc,
                 data = binarydata,
                 studlab = paste(Author),
                 comb.fixed = FALSE,
                 comb.random = TRUE,
                 method.tau = "SJ",
                 hakn = TRUE,
                 prediction = TRUE,
                 incr = 0.1,
                 sm = "RR")
m.bin

其它参数与具有连续结果数据的元分析中使用的参数类似，但有两个例外：

sm：由于想要对二分类数据合并效应，必须选择另一种汇总度量。可以从“OR”（比值比）、“RR”（风险比）或“RD”（风险差）等选项中进行选择。

inc：定义是否以及如何执行连续性校正。如果数据中的一个单元格为零（例如干预组中没有人死亡），则需要进行这样的校正。在某些情况下，这可能是一种常见的现象，并会扭曲我们对效应量的估计。默认情况下，metabin函数在所有n为0的单元格中添加0.5值(Gart & Zweifel, 1967)。可以使用incr-参数（例如incr=0.1）更改此值。如果试验组总n非常不均匀，也可以使用干预组连续性校正(J. Sweeting, J. Sutton, & C. Lambert,2004)。这通过使用incr="TACC"来实现。

输出结果：

L’Abbé Plots

所谓的L’Abbé图(L’Abbé,Detsky, & O’Rourke, 1987)是基于事件率可视化数据的好方法。在L’Abbé图中，研究干预组的事件率与对照组的事件率作图，研究的N由图中气泡的大小表示。尽管它很简单，但这个图允许我们检查元分析的三个重要方面：元分析的总体趋势、效应量的异质性和事件率的异质性。

在metabin函数中能方便的生成L’Abbé图结果，使用meta包中labbe.metabin函数。对于本例，将使用先前使用metabin函数生成的m.bin输出。

4.3.2 Incidence rates

为了使用发病率数据进行元分析，所谓的人员时间（person-time）数据需要手动收集或计算。计算人员时间数据基本上需要的是事件的数量和它们发生的时间范围。使用IRR.data数据集，其中的人员时间分别为实验组和对照组存储在time.e和time.c列。event.e和event.c列包含两个组的事件数（本例中是抑郁发作）。

代码：

metainc(event.e, 
        time.e, 
        event.c, 
        time.c, 
        studlab = paste(Author), 
        data = IRR.data, 
        sm = "IRR",
        method.tau = "DL",
        comb.random = TRUE,
        comb.fixed = FALSE,
        hakn = TRUE)

输出结果：

4.3.3Pre-calculated effect sizes

我们还可以使用metagen函数合并基于二分类结果数据的预先计算的效应量。然而，需要注意的是二分类结果的效应量通常在合并效应量之前先进行对数转换。这在metabin和metainc函数中自动完成的，但是使用metagen函数，就必须自己完成这一步。下面使用bin.metagen数据集：

该数据集包含了7项研究的风险比以及上下置信区间。为了生成metagen函数所需的输入，首先对所有效应量数据进行对数转换。

bin.metagen$RR <- log(bin.metagen$RR)
bin.metagen$lower <- log(bin.metagen$lower)
bin.metagen$upper <- log(bin.metagen$upper)

计算标准误差seTE：

bin.metagen$seTE <- (bin.metagen$upper - bin.metagen$lower)/3.92

现在我们已经拥有了使用metagen函数处理效应量数据所需的一切。为了将效应量重新转换为原尺寸，在函数中指定sm="RR"：

metagen(RR, 
        seTE, 
        studlab = Author,
        method.tau = "SJ",
        sm = "RR",
        data = bin.metagen)

输出结果：

 4.4 Correlations

对相关性进行元分析与之前描述的方法并没有太大的不同。在合并相关性时，建议使用Fisher’s Z转换来获得每个研究的准确权重。meta包中包括一个用于相关性元分析的附加功能，即metacor函数，它为我们完成了大部分计算工作。

metacor函数的参数与前面描述的metagen和metacont函数基本相同（见4.1和4.2章），其背后的统计数据也是如此。使用cordata数据集：

这个数据集可以使用metacor函数计算元分析。与前几章一样，将使用随机效应模型和Sidik-Jonkman估计来评估研究间的异质性τ²。此外，将sm设置为“ZCOR”，以使用相关系数Z转换进行元分析。其余的语法保持相同。

代码：

m.cor <- metacor(cor, 
                 n, 
                 data = cordata,
                 studlab = cordata$Author,
                 sm = "ZCOR",
                 method.tau = "SJ")

 输出结果，它具有与metagen、metabin和metainc函数相同的结构。

从输出中可以看出，这个数据集中的合并相关性是r = 0.35（对于随机效应模型），结果是显著的(p < 0.0001)。该分析的异质性I²很大支持使用随机效应模型。

本章教学视频：

笔记整理：

杨之绪、王丹、张宇航、聂佩欣、赵加伟、王航航、马达、于艺凝

排版：

陈阳

往期精选

R元分析专题介绍：本活动是由荷兰心理统计联盟组织的第三次线上同辈互助学习活动，小组成员每周学习一章，并整理笔记录屏内容共享，从R基础到元分析效应量计算，作图，元分析高阶，如多水平meta, MetaSEM, network Meta-Analysis 。学习资料为《Doing Meta-Analsis in R 》这本书。笔记为小组成员整理（付费与否全由小组成员决定）

欢迎关注我们

一个专注于管理学及心理学学术研究的平台